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Abstract 

This study presents a new triangular finite element 
for modeling thick sandwich panels, subjected to ther- 
mo-mechanical loading, based on a {3,2}-order single- 
layer plate theory. A hybrid energy functional is em- 
ployed in the derivation of the element because of a C 1 * 
interelement continuity requirement. The single-layer 
theory is based on five weighted-average field variables 
arising from the cubic and quadratic representations of 
the in-plane and transverse displacement fields, 
respectively. The variations of temperature and distrib- 
uted loading acting on the top and bottom surfaces are 
non-uniform. The temperature varies linearly through 
the thickness. 

Introduction 

The analysis of sandwich constructions with aniso- 
tropic composite facings for thermal stresses is of pri- 
mary importance for the efficient design of high-speed 
aircraft and missiles. Subjected to hot environments, an 
analysis of such complex vehicles made of sandwich 
constructions can be performed by employing standard 
three-dimensional finite elements. However, this meth- 
od is computationally demanding and very costly. An 
alternative to using standard elements is the employ- 
ment of finite elements based on either a single-layer or 
discrete-layer plate theory in order to eliminate the dis- 
cretization in the transverse direction. Concerning me- 
chanical loading only, Das et al. 1 provided an extensive 
review of the previous studies and introduced a new 
triangular finite element based on a single-layer theory. 

A majority of the previous work concerning only 
thermal loading utilized the first-order shear deforma- 
tion theory (FSDT). Extensive reviews of this theory, in 
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the context of linear and non-linear finite element anal- 
ysis, can be found in Refs. [2]-[5] . Although the FSDT 
provides acceptable mid-plane displacements, as well 
as in-plane stress components, it fails to capture the 
accurate behavior of transverse stresses, especially 
when analyzing sandwich constructions, which may 
have extremely different mechanical and thermal prop- 
erties. 

In order to improve the FSDT, Rolfes et al. 6 devel- 
oped a post-processing technique for computation of 
the transverse stresses by integrating the equilibrium 
equations in conjunction with a predictor-corrector 
approach when analyzing laminates. Park and Kim 7 
also introduced an approach to improve the FSDT by a 
two-phase analysis process. The first phase involves 
computation of the displacement components by inte- 
grating the transverse strains from the finite element 
analysis of FSDT. The second phase involves updating 
the force vector based on the non-linear part of the 
displacement components for further finite element 
analysis. Thankam et al. s developed a material finite 
element for which the displacement approximation in- 
cludes material properties, as well as the element coor- 
dinates. The expressions for the in-plane displacements 
are obtained from modified equilibrium equations. 

Kant and Khare 9 employed a {3,0}-order plate 
theory with cubic variation for in-plane displacements 
and a uniform transverse displacement. This theory 
does not require shear correction factors, and the 
resulting element utilizes a C° isoparametric interpola- 
tion function. However, as discussed by Naganarayan et 
al., 10 this type of element is prone to error in the case of 
thermo-elastic loading. The error is traced to the non- 
participation and/or inconsistent participation of the 
higher-order terms of certain components of the strain 
fields in the total potential. In the case of thick lami- 
nates and sandwich constructions, the major shortcom- 
ing of this theory is the assumption of uniform trans- 
verse deformation. Employing the {1,2} single-layer 
theory and exact solution, Tessler et al." showed that 
the transverse deformation is not uniform. Although the 
single-layer plate theory captures the non-uniform 
transverse deformation, it fails to capture the discontin- 
uous slope of the in-plane and transverse displacement 
components in the thickness direction. 
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The zig-zag theories and layer-wise theories elimi- 
nate this shortcoming of the single-layer theory. These 
theories enforce the continuity of in-plane displace- 
ments and transverse shear stresses at each layer inter- 
face, as well as the zero transverse shear stress at the 
top and bottom surfaces. Ali et al., 12 Kapuria et al., 13 
and Xiaoping and Liangxin 14 extended such theories to 
investigate the response of laminates and sandwich 
beams under thermal loading. Although layer-wise 
theories are suitable for the development of new finite 
elements, they become computationally expensive for 
sandwich structures containing multi-layered face 
sheets. 

In order to take advantage of the computational effi- 
ciency of the single-layer theory, this study presents a 
new triangular sandwich finite element for thermo- 
elastic analysis based on a hybrid-stress finite element 
formulation while employing the {3,2}-order single- 
layer plate theory developed by Das et al. It captures 
both the transverse stresses and the discontinuous slope 
of the in-plane and transverse displacement components 
in the thickness direction by post-processing proce- 
dures. This element is validated by comparison against 
the exact solution for a simply supported laminate and a 
sandwich construction subjected to non-uniform tem- 
perature distribution. 

Element Description 

The sandwich element is composed of two face 
sheets with a soft core in between, as shown in Fig. 1. 
The mid-plane of the triangular sandwich element coin- 
cides with the (x, y) plane of the element coordinate 
system. The origin of the element coordinate system is 
located at a node, and the x-axis is directed along an 
edge. The vector n , represents the unit normal to the 
k th edge. The element has a uniform thickness of 2h. 
The thicknesses of the core and face sheets are denoted 
by h c and h f , respectively. The face sheets, as well as 
the core, can be composed of homogeneous, elastic, or- 
thotropic material layers. Each ply forming the face 
sheets has elastic moduli E n , E 22 , and E 33 ; shear 
moduli G 12 , G 13 , and G 23 ; Poisson’s ratios v l2 , v 13 , 
and v 23 ; and thermal expansion coefficients «, , a 2 
and a 3 , where subscripts 1 and 2 are the material coor- 
dinates. For an isotropic core material, the Young’s 
modulus, Poisson’s ratio, and thermal expansion coef- 
ficient are defined as E c , v c , and a r , respectively. 

As shown in Fig. 1, the sandwich element is sub- 
jected to arbitrarily distributed loading acting on its 
upper and lower surfaces. The load distributions acting 
in the direction normal to the upper and lower surfaces 
are denoted by q + (x,y) and q~(x,y), respectively. 
The concentrated forces in an arbitrary direction are 
also considered as part of the external loading. The top 




Fig. 1 Geometric description of the thick sandwich 
plate element. 

and bottom surfaces of the element are exposed to the 
temperature variation of T + (x, y) and T (x, y) , re- 
spectively, The temperature varies linearly through the 
thickness. Furthermore, the sandwich element is sub- 
jected to arbitrary kinematic boundary conditions. 

As shown in Fig. 1, this triangular element contains 
13 degrees of freedom at each node. These degrees of 
freedom consist of two in-plane displacements (u, v) , 
two out-of-plane rotations (6, . 6 ) , and three transverse 
nodal displacements ( w, w, , w 2 ) and their derivatives 
{w x , w y , w l x , w l y , w 2x , w 2y ) . The weighted-average 
in-plane displacement components in the x- and y- 
directions are denoted by u and v , respectively. The 
weighted-average transverse displacement is denoted 
by w. The weighted-average bending rotations about the 
negative x- and positive y-axes are denoted by 0 , and 
9 , respectively. Their positive sign convention is 
shown in Fig. 1. The transverse displacements, vv, and 
w 2 (not weighted-averaged), represent the symmetric 
and anti-symmetric expansion modes through the thick- 
ness of the element. 

Higher-Order {3,2} Plate Theory 

The {3,2}-order plate theory implies that the in-plane 
displacement components vary cubically and the trans- 
verse displacement component varies quadratically 
across the thickness of the panel. At any point in the 
panel, the in-plane displacement components in the x - 
and y -directions are represented by u x (x, y, z) and 
u v (x,y,z), respectively, and the transverse displace- 
ment component by u z (x, y, z) . In accordance with the 
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{3,2} plate theory, the displacement components of the 
sandwich panel are defined in the form 


u x (x, y, z) = u 0 (x, y) + u ] (x,y)£ 

+ u 2 (x,y)C 2 +u i (x,y)£ 3 


u y (x, y, z) = v 0 (x, y) + v, (x, y)£ 

+ v 2 (x,y)C 2 +v 3 (x,y)£ 3 

u z (x, y, z ) = w(x, y) + w, (x, y% 

+ w 2 (x,y)(c 2 -1/5) 


(la) 


(lb) 


(lc) 


where C, -z/h is the normalized thickness and varies 
in the range -1 <C<1. 

As introduced by Tessler, 15 the weighted-average 
quantities through the thickness of the panel are defined 
as 

1 h 

(u(x,y),v(x,y)) = — \ (u x (x,y,z),u y (x,y,zj)dz (2a) 

-h 

[d x {x,y),d y {x,y)) 

h 

■ J ( m v (x, y, z), u x (x, y, z)jzdz 


J5_ 

2h~ 3 


(2b) 


Mx, y) = -j- f u : (x, y, z)( 1 - C" )dz 

4 h J 


(2c) 


Enforcing the definition for weighted-average kine- 
matic variables results in the expressions for u i and v t 
(i = 2,3) in terms of the weighted-average kinematic 
variables, u, v, 0 X , and 0 V , and the unknown coeffi- 
cients, h, and v, (i = 0,1). The unknown coefficients, 
u, and v, (i = 0,1), are determined by imposing the 
condition of zero transverse shear strains on the top and 
bottom surfaces of the sandwich panel. 

After determining the coefficients it, and v, (z = 0,1, 
2,3), the displacement components, u x and u , are ex- 
pressed in terms of the weighted-average quantities 
(u,v,w,9 x ,9 ) as 


u x (x, y, z) = u(x, y) + h^9 y (x, y) 
\hw lx (x,y) 


6 2 


+ h 


f r /-3V 


V 5 3 

+ w 2 A x >y) 


5 ( 3a ) 

-(0 y (x,y) + w x (x,y)) 


u y (x,y, z) = v(x, y) + h^9 x (x, y) 
\hw Xy (x,y) 


+ 


(1 


6 2 


J 

r 

+ h 

\ 


5 3 


■5 (3b) 

-(e x (x,y) + w y (x,y)) 


+w 2 , y (x,y) 


Substituting for the displacement components from 
Eqs. (3) and (lc) in the strain-displacement relations of 
the elasticity theory, the strain components can be ex- 
pressed in the form 


e a0 = m e a0 +hC m K afi 


( 1 


+ 


I_£l 

6 2 


<»c 


J cc/} 


(a,P = x,y) (4a) 


V” ~ J 
f 

+ h 




V 5 3 y 


0)* 


^ap 


e s = m e a + 2 hC (U, /c 


(ry;,rJ = -a-C 2 )rryJ 0) r x P 


(4b) 


in which the resultant strains and curvatures are defined 
as 


(0) ^ 

= U , 

(0 V = 

= V , 


,x ‘ 

yy 

,y 9 

(0) £ 

zz 

II 

0 

= U y + V x 

II 

2s- 

a 

C*T 

3 

Si 

II 

a 

11 

(0) /C = 

-0 , 

<v 

= 9 

XX 

y,x ’ 

yy 

x,y 

(0)„ - 

K :z ~ 


(»)„ 

= 9 +9 

x,x y. 


(5a) 


- 


= -r(9 )XX + w xx ) + M> 2x 


a) 


K yy = ^ e *,y + w J + w Xyy 
m K xy =^(0 x , x +9 yfr + 2W' Xy )+2w 2 ' Xy 


(O) r y , = w y +0 x , ( °V,z = w,+(9 v 


(5c) 

(5d) 

(5e) 

(5f) 

(5g) 


Unlike the conventional form of the relationship 
between the stress, cr , and the strain, £•.. , components, 
an alternative “mixed” form of the constitutive equa- 
tions is employed as 
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(k) 


(6) 


in which the superscript ( k ) indicates that the quantities 
are obtained from the stress-strain relations in the A' th 
layer, and strain and stress quantities without the super- 
script on the right-hand side of Eq. (6) represent 
through-the-thickness behavior. The components of the 
material property matrix are defined as 


r'W) _ _ n(k)s-i(k) 

ij ij iV /3 3 


W = with sg= — 

Co-, 


(i,j =1, 2, 6) 
(i =1, 2, 6) 


(7a) 

(7b) 


ous through-the-thickness distribution of the transverse 
normal stress, a cubic polynomial is assumed for a zz as 


a zz (x,y,z) = <r z _ 0 (x, y) + </>(£) a . z , (x, y) (11) 

where (/>(%) = £ /3 and <r zz0 and cr z .j are unknown 

functions of x and y. The unknown functions of the 
transverse stress component, <r zz0 and cr /y] , are deter- 
mined such that the transverse strain component £ (k) 
becomes the least-squares equivalent of the kinematic 
definition of the transverse normal strain component 
(i.e., £ zz =u zz ). This is achieved by minimizing the 
total error between £ {k) and s zz through the thickness 
of the panel, i.e., 


n 

J(< 


.(*) 



( 12 ) 


Minimization of Eq. (12) with respect to <r ;z0 and <r zzl 
results in the expression for <t zz0 and o z , , in the form 



(13) 


in which the known vectors q' a and q" (a = s,/c) are 
defined as 


The temperature change from a reference state, AT , is 
assumed to vary linearly through the thickness and it 
has a non-uniform variation defined by T(x,y) such 
that 

AT = (f l + f^)T{x,y) (8) 


ql = Wla 

q :={<?: 


ha 

a 

h a 


C Jla ) 

II 1 

C h a I 


(a - £,k) (14) 


with their explicit expressions given in Appendix A. 
The resultant strain and curvature vectors are defined 
by 


where 7J and T 2 can be expressed in terms of the 
surface temperatures at the top and bottom surfaces, 
i.e., 


(9) 


The transverse normal strain in the k th layer can then be 
expressed as 


" T*{x,y) + T~ Oyy) 

V TiX^) , 

' T + (x,y)-T~(x, y)'' 
, T{x,y) 


- R^Cyy - + Sg<T„ i 

+ ll0) 

in which the continuous through-the-thickness varia- 
tions of s aa and y ap {a, ft = x,y) are known explic- 
itly, but not that of a, z . In order to achieve a continu- 


«= {'»«„. ' V "V'V "V„) (15a) 

K r = {'%,,<%, ‘V,, %.%. "A) (isb) 


Governing Equations 

The governing equations concerning the equilibrium 
equations and continuity of interelement displacements 
along the element edges are derived based on the prin- 
ciple of minimum potential energy. The unknown trac- 
tion components, A x , A y , and A, , acting along the 
element edge ensure the continuity of interelement dis- 
placements. 

For a sandwich element subjected to normal tractions 
on the top and bottom surfaces, c i ( x ,y>h) and 
9 (x, y, -h ) , respectively, and the unknown interele- 
ment tractions, A, , A y , and A. , along the edges, the 
total potential energy, n s , of the element “e” is 
expressed as 
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n = j j (cr^^ + a (k) £ yy +a zz £ ik) +a (k) y J , 

A e -h 

+ cr i- V,, + c^r xy jdzdA- j(q\(x,y,h))dA 

A e 

+| (<1 u z(x, y,—h)jdA (16) 

A e 

h 

- J J [A (W, - ) + A, (»y ~ Uya ) 

r e -h 

+X z {u z -u za )]dzdT 

where A e represents the area of the element mid-sur- 
face and T e the boundary of the element. The displace- 
ment components, u xa , u va , and u :a , belong to the 
adjacent elements denoted by a. In accordance with the 
{3,2} plate theory, the weighted-average quantities 
associated with the adjacent element are represented by 
u a > , w a , d xa , and Q ya . 

The governing equations are obtained by enforcing 
the first variation of the total potential energy with 
respect to the element unknowns (i.e., the displacement 
components, u x , u , and u. , and the interelement 
tractions, A x , A , and A . ) to vanish. This requirement 
results in 

h 

SU = f r ( a {k) 8s + a (k) 8s +a Ss (k) + a (i) Sy 

e J I ( it xx yy yy zz zz yz / yz 

A e -h 

+<?* ) 8r X z + + ^Sy x: + < , d'y K . ) dz dA 

- 1 (<? + Su z (x,y,h)^dA 

A e 

+ 1 (q8u z (x,y,-h)'jdA (17) 

A e 

- J j \SI X (u x - u xa ) + SA y (u y - u ya ) 

r e -h 

+ SA_ (u z - u za )] dz dV 

h 

- J J + A r Su v + KSu z ]dzdT = 0 

r e -h 


+ m QJ m rA dA 

+ 1 - q )5w + (q - q~)5w t H — (q* - q~)Sw, | 

+ \{s m T{u-u a ) + d ,0 ’M (<? - e m ) + S m Th(w t t 


W ) 


+ 5 {) M 


5 W 5 

— (0 +w ) + w - —(6 +w ) + w 
4 J \ 4 


+ 5 m T (v-vJ + 8 m M x (0 - 0 J + 5 '”«(«; y - w, ,J 
+ S 1 1 M (0 + w y ) + w l y j - ( 0 „ + w v ) + w 2a y 

+S ""Q^ (w—w a ) + d ,n Q (w t - w u ) + d''Q n f w 2 — w 2a ) 


+ J J 1 "M <5^— (0 + w ) + w, 

+ ia 'fSv+ '°'M 39 + "'f hSw, 


dr 


+ m M <5 —(0 +w ) + w 
4 


+ QJw+ QJw, + QJwAdT = 0 


(18) 


where the left superscripts (0) and (1) denote the con- 
stant- and higher-order modes of deformations in the 
transverse direction, respectively. The expressions for 
the interelement stress and moment resultants, U] T a , 
u ] Q nz , and w M a ( a = x,y; i = 0,1,2) , are defined as 


( ,<'>f,<'>M x )=\ A x ( P^P^P^P^dz (19a) 


( i0) f y , W M y , m f y , m M y ) = J X y (P n J\J\.P,)dz (19b) 


( m ( L . 0 , a -. P ) a ,) = 


X z \l,^--\dz (19c) 


in which 


Substituting for the displacement and strain components 
from Eqs. (1), (3), (4), and (10) in conjunction with 
Eqs. (11) and (12) and integrating through the thickness 
in Eq. (17) result in 


c-t- r ff<°>Ar e (°) , (0) » r £.(0) , (0), r ^ (0) (0) , T ~ (0) 

o II = < No £ + No £ + No £ + No y 

e J \ xx XX yy yy zz zz xy • xy 


, 0 ),r o(l) , (l) A r £.( 1 ) , (l) A r £.( 1 ) 

+ No £ + No £ + N o y 

xx xx yy yy xy • xy 

, <°> a ^ C*(°) i <°)»^ c*(°) , ( 0 ),, £.( 0 ) C°) » ^ £>( 0 ) 

+ MO K + M O K + Mo K + Mo K 

xx xx yy yy zz zz xy : 


(continued in next column) 


P 0 =l, P \=h$, P 2 = 

E 3^ 




P 3 =h 


5 3 


(19d) 


The resultant stresses , U) N ap , (<>) Q VZ , and m Q xz , and 
the resultant moments, (a,/3 = x,y\ f = 0,l), 

can be expressed in terms of the strain and curvature 
components by integrating the strain energy density of 
the sandwich panel [i.e., the first integrand in Eq. (16)] 
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through the thickness after substituting from Eqs. (4), 
( 6 ), (10), and (11) in conjunction with Eq. (13) as 




"A 

B 

0 
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N* ' 

M 

.= 

B r 

D 

0 

• 

K 



M* 

k 


0 

0 

G 


Y. 


0 


T{x,y) (20) 


where 


-( (1) M 

+ m M ) 

|+ m M 

4 V yy>y 

xy,x J 

' yy,y 


+ i0) K y , x 

i 

■P 2 

II 

o 

~( m M 

+ n> M 

)+ {0) M 

4 \ 

xy,y , 

1 XX, X 


+ 

o 

II 

Ot 

o' 
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(24d) 


(24e) 


N r = { (0) N XX , m N„, ( 0 ) N z: , m N xy , W N XX , ( X, ( X} 

( 21 a) 

(21b) 

N * = { m KrN; y , (0 ) k, w K’ m K> m K> m K} T 

(21c) 

(2 Id) 


Q T ={ (0 %rQ xz 


( 21 e) 


66={ 

o' 

(0) (0) 

XX * 

,i y ’ 1 

ii 

( (C V 

, ,o v , ,c 


1 XX 

’ yy ’ 

II 

t-, 


,‘"V4 


xy 7 xx 7 yy 7 xy 


( 21 h) 


The explicit forms of A, B, D, N , and M are given 
in Appendix B. This constitutive relation can be also 
expressed in terms of the compliance matrix C in the 
form 


(0) W / \ 

—^+h( m N XXtXX + m N myy +2 m N xy>xy ) = q 2 (24f) 

M ° +( W M + m M +2 W M ) = -o. (24g) 

^2 t xx,xx yy.yy ir * xy.xy ) 5 Vl 

where 


< 7 i =c l* ~ C 1 and q 2 = c l + +C 1 

The interelement stress and moment resultants, u) T a , 
( "Q nz , and U) M a (a = x,y; i = 0,1,2) , along the 
element boundaries are obtained as 

'“t-'X-A + 'X", (25a) 

m f y = (0) N yy n y + m N xy n x (25b) 

( 0 ) 4_- = w Q a n x + w Q r n y 

(25c) 

--( n + m M n j 

4 V 1 yy>y n y ~ 1 xy,y n x) 


E = C(S + S,) (22) 

where 


E r ={6, K,y}, S r ={N,M,Q} 
Sf =T(x,j){N*,M*,0} 


(23a, b,c) 


Substituting for the resultant strains and curvatures 
from Eq. (5), integrating the area integrals of Eq. (18) 
by parts, and applying the Gauss theorem to appropriate 
terms result in the governing equations for the sandwich 
element. The equations of equilibrium are obtained as 


m M y = w M„,n y + {(>) M xy n x 

(25d) 

m M x = m M xx n x + m M xy n y 

(25e) 

] Q„=-h( W N^ x n+ m N my n y 

+ m N xyiX n y + m N xyty n x ) 

(25f) 

m t = m N xx n x + m N xy n y 

(25g) 

II 

p 

s 

V 

+ 

S 

(25h) 


m N + m N = o 

xx, x xy,y 

(24a) 

af 

si 

+ 

s4 

I? 

1 

II 

(25i) 

(0) N + im N = 0 

yy,y xy,x 

(24b) 

+ m M my n y + m M xyy n x ) 

d + m M +2 m M ) 


( 1 ) M, = a) M xx n x + 

(25j) 

xx, xx yy,yy xy,xy j 

- W Q - (0) Q = q. 

2oxz,x ^yz,y 1 1 

(24c) 

W M y = m M yy n y + m M xy n x 

(25k) 
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The interelement kinematic continuity conditions are 
obtained as 

u=u a , v = v a , w = w a ; 9 X = 6 m , 9 y = 9 ya (26a) 

W 1 = W l„ > W l,z = W la,x ’ W l,y = W la,y ’ W 2 = W 2a ( 26b ) 


displacements, respectively. The compliance matrix, 
C , is defined in Eq. (22). In accordance with the hy- 
brid energy formulation, the resultant stress vector, S , 
must satisfy the equilibrium equations, Eq. (24), iden- 
tically. It is composed of the homogeneous solution, 
S ; , , and the particular solution, S a , in the form 


r 5 . 


5 

+ 

+ 

i ^ 


-(6ya +W a,x) + W 2a,x 

r 5 i 


"5 „ 

4 (f/, +iv,)+ny,. 

— 

-(&«,+ W a,y) + W 2a,y 


(26c) 

(26d) 


These kinematic continuity conditions are imposed not 
only on the weighted-average displacements and slopes 
(u,v,w,9 x ,9 ) , but also on the derivatives of the 
higher-order displacement modes (w l ,w 2 ) in the trans- 
verse direction. Therefore, the finite element implemen- 
tation of the equilibrium equations requires at least C 1 
interelement continuity for the out-of-plane displace- 
ment modes of w, vi , , and w 2 . Because of this require- 
ment, the finite element implementation of the total po- 
tential energy functional in terms of the assumed dis- 
placement field becomes rather difficult. Therefore, this 
triangular sandwich element is developed by utilizing 
the hybrid energy functional in which the stress/ 
moment resultants within the element and displace- 
ments and slopes around the element boundaries are 
assumed independently. 


S = S,+S 0 (28) 

The components of the homogeneous part of the 
resultant stress vector, S ; , , are expressed in the form of 
complex power series as 


4 


{0) N =V (0) N r (0) a 

JV ^ Zj L 'Zl(k) U ^(*) 

k= 0 
4 

(29a) 

(1) N =V (1) N r (1) a 

k= 0 
2 

(29b) 

o' 

■Ai 

z 

o 

ii 

z t! 

o' 

(29c) 

k= 0 

4 

(29d) 

{l) M = V (1) M r (1) 3 

lvl 4n iyi 4rj(k) P fr(k) 

k-0 

2 

(29e) 

(0) M =V (0) M r (0) B 

JV1 zz 1TA zz(/t) Pzz(/t) 

k=0 

(29f) 


Element Development 

The hybrid energy functional formulation overcomes 
the difficulty of the C 1 interelement continuity require- 
ment because the displacements, as well as the slopes, 
are independently assumed only along element bound- 
aries, which can be rendered identical along the com- 
mon boundaries of adjacent elements. However, the 
kinematic compatibility between the displacements and 
slopes along the element boundaries is preserved in 
order to avoid a possible shear-locking phenomenon. 
Also, as part of the hybrid energy functional formula- 
tion, the stress and moment resultants within the ele- 
ment are selected in the form of complex power series 
so as to satisfy the equilibrium equations of Eq. (24). 

The hybrid energy functional for an element, 11 H , is 
defined as 

n„ =-Js r CSt/4+jsfEt/A-jT>„t/r (27) 

2 a , A e r e 


<0, 2fz = Z (0, Qfz(*, (0, P f z(A-) (29g) 

k= 0 

in which the vectors of stress, (<) N* , moment, , 

and shear stress, (0) Q* Z , resultants contain the real and 
imaginary parts of the complex power series terms, and 
the vectors and l0) p iHk) contain the 

real and imaginary parts of the unknown coefficients 
(generalized coordinates) with $,ri = x,y, * = 0,1. The 
series representations are truncated at an appropriate 
order so as to ensure compatibility among the stress, 
moment, and shear stress resultants, when substituted 
into the equilibrium equations, and to ensure definite- 
ness of the stiffness matrix associated with the general- 
ized coordinates. With this representation, the homo- 
geneous equilibrium equations are satisfied, provided 
that certain relations exist among the unknown coeffi- 
cients. These constraint conditions are expressed as 

cb = 0 (30) 


in which the element boundary is denoted by E e and its 
area by A e . The vectors T () and u i; include the 
components of the boundary forces and boundary 


in which c is the constraint matrix and b is the vector 
of unknown coefficients. The constraint equations that 
are required to construct matrix c are given in Appen- 
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dix C. The vectors (<J N^ (t) , (0 M fl7(Jt) , and (0, Q flW , 


and their 

associated 

unknown 

coefficient vectors, 

w a 

k ) ’ 

: °P *,(*>» and 

(0) o 

are expressed in the 

form 




(0 N r 

= ( °M r 

= {Re[z*],. 

- Tm[z* ]} 




{k = 0, 1,2, 3,4) (31a) 

(0 N r 

Ly zz(k) 

= = 
1V1 zz(lr) 

{Re[z 4 ],- 

' hn[z A ]} 




(k = 0,1,2) (31b) 

(°>o r 

= {Re[z 4 ],- 

lm[z*]} 





(k = 0,1,2, 3) (31c) 

<f) a r 

sW) 

= {Re[ (<) « W )]Jm[ (f) a^ (t ,]} 




(k = 0, 1, 2,3,4) (32a) 

<n a r 

U'zzfjt) 

= {Re[ < V-(i ) ], I m[ (0 « 

Z -(*) 


(k = 0,1,2) (32b) 


(k = 0,1, 2,3,4) (32c) 

( °PL W ={Re[ (0 ^ (i) ],Im[< < ^ ) ]} 

(^ = 0,1,2) (3 2d) 

(0 V^)={Re[ ( % m ],Im[ ( V w ]} 

{k = 0,1,2, 3) (32e) 

with ^,rj = x, y and ^ = 0,1; the complex variable, z , 
is defined as z = x + iy , in which / = V-T . The vector 
representations of the stress, bending, and shear stress 
resultants can be cast into matrix form as 

S*=Pb (33) 

in which the matrix P contains the real and imaginary 
parts of the power series terms and the vector b con- 
tains the real and imaginary parts of all the unknown 
generalized coordinates. The explicit expressions for P 
and b are given in Appendix C. Hence, the general 
solution form of the resultant stress vector S becomes 

S = Pb + S D (34) 

As part of the hybrid finite element formulation, the 
displacement field is assumed along the element bound- 
ary. Along the k th edge of the triangular element lying 
between nodes i and j (as shown in Fig. 2), a Carte- 
sian coordinate system is attached to node i and the 
axes normal and parallel to the £ th edge are referred to 
as h and s, respectively. The .s-axis also represents the 
line coordinate along the boundary, and it varies be- 



Fig. 2 Assumed displacement and slope fields along 
the element boundary. 

tween 0 < s < l k , with i k being the length of the k th 
edge. Therefore, the boundary displacements and their 
derivatives, as well as the independent slope compo- 
nents, as defined by the {3,2} -order theory, are expres- 
sed in terms of the nodal degrees of freedom as 

(«<« ) = X iV* (Jf) ) 

m = 1 

(35a) 

2 

( ,.,(*) ..,(*) j-r fcA/'vt ,(*) w (*) ,„(*) \ 

m = 1 

+ H m+2 ( S )(^(l) ’ W U(») > W 2,s(m) ) 

(35b) 

i 

2 

> W l,n ’ W 2,n ) / , m W 1 W .n(m)’ W ].n(m^’ W 2.n(m) ) 

m= 1 

(35c) 

9 (k} = Vy (s)0 <k \ 

s / j m V / s(m ) 

m= 1 

(35d) 

o { n k) =i M ,ss) c, 

m-l 

(35e) 


where ( u (k \u[ k] ) represents the in-plane displacements 
normal and tangent to the k lh edge and ( , w \ k] , w {k) ) 
and ( w [k) , wf k J , ) denote the weighted-average, 
first-, and second-order modes of out-of-plane dis- 
placements and their derivatives, respectively. The 
slopes normal and tangent to the boundary are denoted 
by 9 (k) and 9 (k) , respectively, and their sign conven- 
tions are shown in Fig. 2. The corresponding nodal 
degrees of freedom for these displacement and slope 
components at nodes i and j along the k ,h edge are 
identified by the subscript m=\ and m= 2, respectively. 
Note that the slope component, 9\ l: 1 , has a mid-side 
(third-node) degree of freedom to be eliminated later. 

The functions N.(s), M (s ) , and Hj(s) represent 
the linear, quadratic, and cubic (Hermitian) interpola- 
tion functions, respectively. Also, the out-of-plane dis- 
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placement components, w (k) , w \ k) , and w [ k) , are func- 
tions of the 5-axis only. Therefore, the derivatives of 
( w^ k) ,w[ k) ,w ^ k) ) with respect to the /?-axis are assumed 
independently by using linear shape functions. 

In order to maintain kinematic compatibility, differ- 
ent orders of polynomials are chosen to represent the 
displacement and slope components as given in Eq. 
(35). The degree of freedom at the mid-side node, 0 ^ , 
is eliminated by imposing the constraint of linear 
variation of the boundary expression 


d 2 

ds 2 


^K ) +0 + < ) 


= 0 


(36) 


With this constraint equation, the mid-side degree of 
freedom, 0^1 , can be expressed in terms of 0 {k) . , 


„(*) 


n( 31 ’ . , 

r \Hm) ’ and w 2 ,l ’ wd h ( m = 1, 2 ), which are assigned at 
the ends of the k th edge. 

In compact form, the in-plane displacements, u k) 
and u (k> , the transverse displacements, w (k) , w [ k) , and 
w (k) , their derivatives with respect to the /?- and 5- 
coordinates, and the slopes around the n- and 5-axes, 
0 ik) and 0 (k) , along the k th edge of the element can be 
expressed as 


„'<*) = L' w v' (t) 


where 


(37) 




= { u l k) ’ u l k) ,w {k \w\ k) ,w {k) ,w[ k) ,w[ k J ,w {k l , 


(38a) 


»(*r _]„'(*/ ,,'(*)■ 


>(kY 


,W ' ,W 


r(kf r (k) r i(k? 
" .« 5 ” 1,« ’ w 2 ,n 


,W 2 

'( kf 




W 


'<*) w '(*)' w 'W tfikY 


(38b) 


in which u 


'(*) ^ /,,,(*) ,.,(*) i . q »(*r —saw am i w ' 


’<kf _/..(*) .,(*) 1 '(*/ _(,,,(*) 

l W a(l)’ M a(2)J > W l W (l) ’ W (3)-< ’ 


- 


W 

M*W ( %}, and ^i!a Y ={ W iMD’ W iM2)}’ with (* = 
1,2; a = n,s). The explicit form of the matrix of 
interpolation functions, L' (A| , is given in Appendix C. 
The prime on the superscript denotes definition with 
respect to the line coordinate system, (n,s ) , along the 
boundary of the element. 

In order to utilize this equation in the hybrid energy 
functional, the vector v ,(k> is expressed in terms of the 
element nodal vector v as 


(*) 


= Z ik) \ 


(39) 


where Z tk) is the transformation matrix between the 
nodal degrees of freedom defined in the element. 


(x,y), and line coordinate, ( n,s ) , systems for the k th 
edge, and v T = {v [ n , \ T {1) , v[ 3) } with 

<A-, = {«(*), v (t) ,w {t) ^ 

W 2,x(k) ’ W ,y(k) ’ W l,y(k) ’ W 2, y(k) ’ ^x(k) ’ ^y(k) } 

The form of the transformation matrix, Z a) , is given 
explicitly in Appendix C. The displacement vector 
u' a> becomes 


u' ( *) =L' (i) Z (i) v (41) 

A similar transformation between the displacement vec- 
tor, u' <k> , defined in the line coordinate system and the 
displacement vector, u lA) , defined in the element 
coordinate system is established as 

u'W =y ( V h (42) 

where 


,(*> 




(43) 


and the transformation matrix Y U| is expressed in 
terms of the sine and cosine of the angle between the 
normal (/7-axis) at the k' h edge and the x -axis. The 
explicit expressions for the components of Y <k> are 
given in Appendix C. It must be noted that the sign 
convention of d[ k> is in the negative x-direction of the 
element, whereas the slope 6 x(k] at the nodal point is in 
the positive x-direction of the element coordinate 
system. Also, the transformation matrix is orthonormal, 
i.e., Y (/:) = Y (A) . Hence, the displacement vector, 
u (/ ) , in the element coordinate system is expressed in 
the terms of the nodal unknowns (or degrees of 
freedom) in the element coordinate system as 

u (t) =L W (5)V (44) 

in which L (A, (5) is defined as 


L w (5) = Y wr L' (A) (5)Z (A) (45) 

The boundary displacement vector, u^ 0 , along the k th 
edge of the element is defined as 




= <u ,v' 




< , T(»;* , +< | ) (46) 

It is related to the displacement vector, u (A) , in the 
element coordinate system as 
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u«=B u u ( * , =B u L w ( S )v 


(47) 


in which B a , given explicitly in Appendix C, can 
easily be constructed by comparing the components of 
the vectors u l/f) and uJ, A) . The boundary stress vector 
T f b k) , containing the resultant stresses and moments 
corresponding to the boundary displacement vector 
Uj A) , can be obtained from 





T l k) = B t S (A) (48) 

in which B, is a matrix transforming the resultant 
stress vector defined in the element to the resultant 
stress vector defined along the k' h edge of the element 
and the explicit form of B s is given in Appendix C. 
While invoking the constraint equations, Eq. (30), and 
substituting for the stress vector, the boundary 
displacement and boundary stress vectors in the hybrid 
energy functional result in 

U H =-b r Hb + k r cb + A R>- R> 

H 2 00 (49) 

+ 'Rf b-b r Gv + n 0 

where 


H = | P r CP dA 

A e 

(50a) 

£ J p r b‘ A) b; A) l (A) dr 

k r m 

(50b) 

" Z j s » K’K ’ L<i) dr 

k=1 r (« 

(50c) 

b R a = \s T 0 CVdA 

A e 

(50d) 

b R, = Jsf CP dA 

A e 

(50e) 

-js T 0 CS 0 dA+ Jsf CS, dA 


' A e A e 

(50f) 


+ Jsf CS 0 dA 


A e 


The unknown vector, l , contains the Lagrange multi- 
pliers used for enforcing the constraint equations, Eq. 
(30). In matrix form, the hybrid energy functional, L1 H , 
can be rewritten as 


n H =ib r Hb + ‘R>- R> 

2 ° (51) 

+ fc Rf b+b r Gv + n, 

where 


In accordance with the concept of energy minimiza- 
tion, the first variation of the hybrid energy functional 
with respect to the unknown vector b of generalized 
coordinates yields 

<?b r (Hb + A R o + A R, - Gv) = 0 (52a) 

or 

b =H“‘ |Gy- *R 0 - 4 R, j (52b) 

With this explicit solution form, the hybrid energy 
functional becomes 

n„=-Vkv+f:v+f, r v+n, (53) 

in which the stiffness matrix k and the resultant force 
(load) vector f are defined as 

k = G r H 1 G , ff = h R T 0 H 1 G- r R^ (54a) 
and 

f, r = *Rf H- 1 G (54b) 

Associated with applied tractions, q + (x,y) and 

q (x,y) , on the upper and lower surfaces of the ele- 
ment, the load vector requires the explicit form of the 
particular solution, S 0 , to the equilibrium equations. 
This requirement is circumvented by computing the 
potential energy of external forces from 

K vj [q + (x,y)uJx,y, + h/2) 

A (55) 

+ q~ (x, y)u z {x,y,-h/2)jdA 

In order to evaluate the integral in Eq. (55), the 
weighted-average variables arising in the expression of 
the transverse displacement are approximated in terms 
of linear interpolation functions (only for the 
computation of the load vector), N m (x, y) , as 

(w,w u w 2 ) = ( 56 ) 

m = 1 

With this approximation, the load vector is recast as 
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f. r = {CM) (57) 

where i T ok = {0, 0, , f Wi J Wiit> , 0, 0, 0, 0, 0, 0, 0, 0} , with 

L m = j N k [q + (x, y)-q-(x,y)\dA (58a) 

A e 

f nm = J N k [_q + (x, y) + q~(x, y)\dA (58b) 

A e 

A (ti =T \N k \_q + (x,y)-q-(x,y)\dA (58c) 

^ A e 

Finally, the element equilibrium equation is obtained by 
requiring the first variation of the hybrid energy 
functional to vanish 

( ?n„=Jv r (kv-f ) -f,) = 0 (59) 

For arbitrary variation of 8x , the element equilibrium 
equations become 

kv=f,+f, (60) 

Assembly of the element stiffness matrices and force 
vectors leads to the structural equations of equilibrium 
as 

K V = F + F. (61) 

in which K , V , and F represent the assembly of the 
element stiffness matrix and nodal and force vectors. 
The vector F arises from the presence of concentrated 
forces. 

Although the transverse shear stresses can readily be 
obtained from the stress-strain relations given in Eq. 
(6), the presence of discrete material properties in the 
thickness direction causes discontinuous variation of 
transverse shear stresses. Flowever, the true variation of 
the transverse shear stresses must be continuous 
through the thickness of the sandwich panel. The accu- 
rate computation of in-plane stresses and their first 
derivatives suggests that continuous transverse shear 
stresses, a xz and a , can be obtained by integrating 
the in-plane equilibrium equations in the thickness 
direction. Similarly, the transverse normal stress com- 
ponent, <7„ , is computed by integrating the transverse 
equilibrium equation in the thickness direction. 
Sandwich panels exhibit nonlinear deformation, with a 
discontinuous slope in the thickness direction due the 
severe changes in mechanical and thermal properties. 
Flowever, in single-layer theories, through-the-thick- 
ness variation of the displacement components is 
approximated by polynomials. Such kinematic repre- 
sentation does not necessarily capture the accurate 
behavior of displacement components in the thickness 
direction because of the presence of discontinuous 


slopes at the interface of the core and face sheet. This 
inaccuracy is avoided by computing the displacement 
components from the integration of the strain compo- 
nents through the thickness. Details of the integration 
procedure for the computation of transverse shear, nor- 
mal stress, and displacement components in the thick- 
ness direction are given in Appendix D. 

Numerical Results 

In order to demonstrate the capability of the present 
element, an unsymmetric laminate and a sandwich 
panel previously considered by Tessler et al. 11 are 
studied. The validity of the present analysis results is 
established by comparison against the three-dimen- 
sional exact solutions presented by Tessler et al. The 
exact solutions were obtained using the PAND3D 
analysis program developed by Burton and Noor. 16 

The square geometry of the unsymmetric laminate 
and sandwich panel are described by 2a = 2b = 20/? , 
with h = 0.5 in. The laminate consists of carbon/epoxy 
plies of equal thickness and a stacking sequence of 
[0 4 /90 4 /0 4 /90 4 ]. In the case of the sandwich panel, 
the face sheets at the bottom and top of the core are of 
equal thickness, h f = 0.2 h ; hence, the core thickness, 
h c , is 80% of the total panel thickness. The face sheet 
material is composed of carbon/epoxy plies, with 
stacking sequence [0 4 / 90, / 0 4 / 90, / 0 4 ] . Each carbon/ 
epoxy ply has Young’s moduli of E n = 22.9xl0 6 psi 
and E 22 = E 2i =1.39xl0 6 psi , shear moduli of G u = 
G u =0.864 xlO 6 psi and G 23 = 0.368 xlO 6 psi , Pois- 
son’s ratios of v n =v 13 =0.32 and v 23 = 0.48, and 
thermal expansion coefficients of or, = 0.01 x 10 6 / °F 
and a, = a 3 = 11.5xl0 -6 / °F . The core is made of soft 
PVC (foam) with isotropic material properties specified 
as E c = 15.x 10 3 psi , v.=0.3,and a c = 15.64 x 10 -6 / 
°F. The shear correction factors that appear in matrix 
G of Eq. (B9) are c v: =c xz = 0.8939 for the unsymmetric 
laminate and c =0.6161 and c E = 0.5536 for the 
sandwich panel. 

In accordance with the {3,2}-order plate theory, the 
boundary conditions are enforced as 

u = w = w 1 = w 2 = w x = w t _ x = w 2x = 0 V = 0 
for y = 0 and y = 2b , and 

v = w=w l =w 2 =w=w l =w 2 =0 JC = 0 
for x = 0 and x = 2 a 

As shown in Fig. 3, the non-uniform temperature dis- 
tributions on the top and bottom surface are, respec- 
tively, specified by T*{x,y) = 300T(v, y)°F and 
T~(x,y) = (>0T(x,y) o ¥ , where T(x,y) = sin(^v/2a) 
xcos(^'>’/2Z)) . A total of 2500 elements are employed 
in the finite element discretization, leading to 16250 
degrees of freedom. 
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h y 




(b) 


Fig. 3 Non-uniform surface temperature applied to 
(a) the |0 4 /90 2 /0 4 /90 2 /0 4 ] laminate and (b) the 
sandwich panel. 

The deformed configurations of the laminate and the 
sandwich panel are shown in Fig. 4. The laminate ex- 
periences significant bending and the upper surface 
bulges out more than the lower surface. In the case of 
the sandwich panel, there is no significant bending of 
the lower face sheet; however, the upper face sheet ex- 
hibits significant bending or it bulges out. The unsym- 
metric stacking sequence of the laminate enables the 
coupling between the in-plane stretching and transverse 
bending due to the non-zero values of the components 


Undeformed configuration 




Fig. 4 Cross-sectional view of the deformed lami- 
nate and the sandwich panel subjected to thermal 
loading. 

of matrix B of Eq. (B3). The variation of u x and u, as 
observed in Fig. 5 shows non-linear behavior, and the 
finite element solution is in close agreement with the 
exact solution. The nonlinear variation of the displace- 
ment components suggests the necessity of higher-order 
theories for thick laminates. Comparisons of the me- 
chanical strains e xx = e xx - a xx AT and s X7 = e :z - 
a,,AT are shown in Fig. 6. The magnitude of the in- 
plane normal strain is comparable to that of the trans- 
verse normal strain due to the steep variation of the 
temperature in the thickness direction. The solutions 
from the present analysis for the in-plane normal stress 
a xx and the transverse shear stress a v . are in close 
agreement with the exact solution, as observed in Fig. 
7. The process of integrating the equilibrium equation 
to obtain the transverse shear stress proves to be effec- 
tive enough, and the solution satisfies the zero shear 
stress conditions at the top and bottom surfaces. Con- 
sidering the complex, meandering nature of the trans- 
verse normal stress a shown in Fig. 7c, the finite 

element results are in acceptable agreement. In the case 
of lower-order theories, such as the {1,2} theory, the 
variation of a,, becomes less accurate for laminates 
with low span-to-thickness ratios, as shown by Tessler 
etal. 11 

In the case of sandwich panels, the presence of a soft 
core in between stiff face sheets gives rise to a highly 
nonlinear behavior of the displacements. As observed 
from Fig. 8, apart from being nonlinear, the in-plane 
displacement u x has a slope discontinuity at the inter- 
face of the core and face sheets. The close proximity of 
the present solution to the exact solution is due to the 
utilization of the process of displacement recovery from 
integration of strain components. The variation of u, is 
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Fig.5 Through-the-thickness in-plane and transverse displacements, u x and u , , for a |0 4 /90 2 /0 4 /90 2 /0 4 ] 
laminate. 




Fig. 6 Through-the-thickness mechanical in-plane normal and mechanical transverse normal strains, s™ 
and £•" , for a |0 4 /90 2 /0 4 /90 2 /0 4 ) laminate. 
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Transverse normal stress o..(a,b,z) psi 

Fig. 7 Through-the-thickness in-plane normal, transverse shear, and transverse normal stresses, <r M , a xz , 
and cr„, for a [0 4 /90 2 /0 4 /90 2 /0 4 | laminate. 




Fig. 8 Through-the-thickness in-plane and transverse displacements, u x and u , , for a sandwich panel. 
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fairly accurate, as illustrated in Fig. 8b. The results of 
the mechanical strains, sf and sf , shown in Fig. 9, 
reveal that even in this case the transverse normal stain 
is comparable to the normal in-plane stain. As depicted 
in Fig. 10, <T X1 and are in close agreement with the 
exact solution. The transverse shear stress follows a 
zig-zag path in the face sheets and the finite element 
solution is able to capture this pattern very accurately, 
apart from satisfying the zero shear stress criteria on the 
top and bottom surfaces. The transverse normal stress 
( 7 ZZ has a parabolic nature and, even though the initial 
assumption of a, z in Eq. (11) is devoid of any quad- 
ratic terms, the equilibrium equation integration process 
recovers the correct solution although the magnitude of 
a zz in the core area is slightly underestimated in the 
present analysis. 

Conclusions 

A new sandwich element based on the {3,2 (-order 
plate theory in conjunction with the hybrid energy func- 
tional has been developed and validated for thermo- 
elastic analysis of composite structures with sandwich 
constructions. The validation was established by com- 
parison against the exact solution for a simply sup- 
ported, unsymmetric, graphite/epoxy laminate and a 
sandwich panel made of a soft PVC core with graphite/ 
epoxy face sheets subjected to a non-uniform tempera- 
ture field. This new computationally efficient element 
captures the correct variation of deformations in the 
thickness direction, as they are evaluated from a pos- 
teriori calculations of equilibrium equations. 
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“ , for a sandwich panel. 
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Fig. 10 Through-the-thickness in-plane normal, transverse shear, and transverse normal stresses, cr a , cr v . 
and cr„ , for a sandwich panel. 
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Appendix A 

The vectors q' s , q", q ' . and q '' in Eq. (1 1) are obtained from the matrix operation 


< qZ 


A, 

Su 

-1 

A 

A 

nQ 

05 ^ 

1 


_s n 

■V 


A 

A 


(Al) 


where S,. (i,j= 1,2) , t’ a , and t" (a = s,k) are given as 

K h k 

[S n S u ,S 22 ] = Y,\[ X ’M 2 ] S K f dz k 


*=1 h t _ 


t a={C} atld *«={£ 


(A2) 

(A: = 1,7; « = £,/c) (A3) 


in which 


tje = Z j < P J s^j 3 S 33 > dz k , t’ jK =Z j (p jK R$S% l dz k 


k=\ h 


*=1 K_ 


A = Z j <pJR%S$dz k , t; =X f VjJRfSgdzt 

*=> h k _ i 


*=1 /i. 


(j = 1,7; no summation on /) (A4) 


with 


and 


1 

/v /3 j \ JV 13 ’^23 > 1 >' fv 63 ’^13 ’ Jv 23 > ''(,3 J 


M = { p oAA^AAA} > {<p ] .} = { p ^ p x^ p iA’ p iA} 

The computation of r t in Eq. (13) is done in the following manner: 


(A5) 

(A6) 


f r >} = 

A, 

A 2 ’ 

"m 

UJ 

Al2 

■v 

W 


where 


/, = Z f (AAAA +Rgsga™+Sg>a? +R£s£a%>)(T 1 + T 2 C)dz k 


* =1 v< 


K “k 

h = Z f K ( 3 AXf +R ( 2 W + A'fAf +^3 , A3 , < , )(r i +T 2 c)(/>dz 


(A7a) 

(A7b) 


*=1 fe 


Appendix B 

The matrices A, D, and B in Eq. (20) are defined as 


A = 


sym. 


A A A4 a 15 A 16 a 17 

A 2 ^23 A 4 A 5 As A 7 

A 3 A 4 As As A 7 

.4 A A A 

^44 ^45 ^46 ^47 

A5 As6 A7 

As A 7 


(Bl) 
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Ai 

D n 

D n 

d 14 

Dn 

A 

A 



£>22 

A 3 

D u 

d 25 

As 

D 




Dn 

Du 

As 

As 

D 





d 44 

As 

As 

D_ 






As 

As 

D 

sym. 





As 

A 

L 







D. 


A, 

Bn 

Bn 

B u 

Bn 

As 

Ay 


A, 

B n 

b 23 

B 2 4 

As 

As 

Ay 


Bn 

Bn 

A 3 

Bn 

Bn 

As 

Ay 

= 

Bn 

B 22 

Bn 

B 44 

As 

As 

Ay 


Bn 

Bn 

Bn 

Bn 

As 

As 

Ay 


Bn 

B 62 

Bn 

A 4 

As 

As 

Ay 


B n 

B 72 

Bn 

A 4 

As 

As 

Bn 


and the matrix G has the form 



(B2) 


(B3) 


(B4) 


By integrating the strain energy density of the sandwich panel [i.e., the first integrand in Eq. (16)], through the 
thickness after substituting from Eqs. (4), (6), (10), and (11) in conjunction with Eq. (13), the components of 
matrices A, B, D, and G are obtained as 


and 


4 =Z f {^ij )( Pis<Pje +-SW4&* 

k=l 4-, 

B j =Z f • 

k=1 4-1 

D ij =Z f {c5*V|r^r+*SVte^r}^* 

* =1 4-i 




c, 


/+3, 7+3 


A- 


(/, / = 1,7 ; no summation on i and /) (B5) 


(*,7 = 1,2) (B6) 


in which CA is defined in matrix form as 



Ai 

^12 

0 

As 

A, 

c 

'-'12 

As 


^22 

0 

As 

Ay 

r 

^22 

As 



0 

0 

0 

0 

0 




As 

As 

c 

y -'26 

As 

sym. 




A, 

r 

M2 

As 






c 

^22 

As 







As 


(B7) 
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and the functions <p. and (p jK ( j = 1,7) are defined in the form 


V'j'=4j t +fa'jA 
¥j K =q'j X +H"j K J 


(./ = U) (B8) 


Note that matrices A, D, and G are symmetric, whereas matrix B becomes nonsymmetric due to the presence of the 

Vj ‘s. 

Although the single-layer theory predicts the correct through-the-thickness variation of transverse shear strains 
for an isotropic panel, it does not consistently predict accurate results for composite sandwich panels having 
significant differences in material properties between the core and the face sheet. The commonly accepted approach 
referred to as the energy (transverse shear strain energy) matching technique is invoked in order to correct for the 
transverse shear strain behavior. This correction is invoked in matrix G as 


G = 


'fin 


c c 

yz xi 


C yz C zfi\2 

dG-,-, 


(B9) 


where c and c x _ are known as the shear correction factors. The shear correction factors improve the global 
transverse displacements of the sandwich panel. The exact transverse shear strain energies are related through the 
shear correction factors as 


U^° a) = -X- U^ correcled) (a = x,y) (BIO) 

C az 

in which c az , with (a = x,y ) , represents the transverse shear correction factors. 

Determination of the transverse shear correction factors for a particular material system requires the exact 
solution of benchmark problems for which the transverse shear strain energy must be independent of the geometry 
and loading. This leads to the shear correction factor being dependent on a particular material system only. 

A square, simply supported, thin panel with a length-to-thickness ratio of 10000 under a double sinusoidal 
distributed loading is considered as a benchmark problem. Based on the analytical solution of the present single- 
layer theory provided by Bamt et al. 17 and the exact solution developed by Anderson et al., 18 the transverse shear 
strain energies u^ uncorrec,ed) and of the panel can be expressed in the form 


2b 2a / 

( uncorrected ) jj ( uncorrected) J* J* ^*^2 I 

0 0 V 

xx_\ 
2 a J 

f . 2 f 

V/7 

^) dxd} ’ 

(Blla) 

2b 2a / 

U ( ‘ xac 0 = J J Sin 2 — 

oo \2a 

V > 

LS//7 

( *y 
V 2b 

j dxdy 

(Bllb) 


where 2 a and 2b are the lengths of the panel along the x - and y -axes and jj ( “ ncorrec,ed ) anc j jjG™*'* are gi ven by 

N z k + i 2 

jj(uncorrec,ed) = £ J Q* ^(uncorrected) ( z )^ z (B 12a) 

k= 1 z i 


u: 


(exact) 


N z k + 1 




-dz 


(B12b) 


where G x =G n and G - G 22 . Therefore, substituting from Eqs. (Bll) and (B12) into Eq. (BIO) results in the 
transverse shear correction factors, c 2 az , being dependent only on the material system in a simplified form 


t r (uncorrected) 

d=- 


t r(exact) 
U Yaz 


(a = x,y) (B13) 
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(B14a) 


The components of the resultant thermal force vector N in Eq. (20) are defined as 


i0> K = “Z I [(?U + Hi){'] + H 2 % 


(*) 

33 


4=1 4-1 


- )(t; + t^)] 

<0>Ar vv = -Z 1 [fe +HL)( r i +fa)s * 

4=1 »*-i 

- («« + <cg> + a%C % ) (t; + t 2 c)] dz t 

m K = -z } [(^. +fau){ r i +< /> r 2) s ^y z k 


4=1 h- 


= ~Z j [(dl+H. L)('i +^ 2 )^ 


(*) 

33 


4=1 4-1 


- K’C + <C<*> + )(7] + T 2 c)\ dz k 


WN >-i | [UL + f^LH'i +H 2 )s\ 

4=1 4-1 


(*) 

33 


+a%ci? +< ) c<* ) )(r 1 +7;^) 


6 2 


fife,. 


( " N l = ~Z j [(?L +HL)(f\ +<f»- 2 )s : 


(k) 

33 


4 = 1 4-1 


6 2 


-(«M +<c4 t) +< ) ^ , )(r 1 +r 2 <r) 

(1) ^ =-Z} [(?7* + H' lK ){r l +</>r 2 )s£ ) 

4=1 4-1 

-(«M +<c<:> +< ) c<* , )(t; + f 2 ^) 

The components of the resultant thermal moment vector M in Eq. (20) are defined as 

W K X = ~Z 1 [( sC + HI ) ('i + H 2 )s : 


dz. 


6 2 


dz. 


?(*) 

J 33 


4=1 4 - 


-^(a£>c?» +<C +< ) c£>)(r 1 +r 2 ^)]<fe t 

<0,M ^ = ~Z 1 [(?** + H'L )('i + h 2 K; 


?(*) 

->33 


4=1 4-1 


- + a%C% + <C<*> ) (t; + T 2 C)] dz k 


(B14b) 


(B14c) 


(B14d) 


(B14e) 


(B14f) 


(B14g) 


(B15a) 


(B15b) 
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K k 

( °K = -X j [{q’ e + HL )('! + <t>r 2 ) Sg> ]dz 


*=i in 


(B15c) 


;<*> 

J 33 


(0>A C = “X j [(^ + ML )( r i + M M: 

k=1 h - 1 

- kT (a™C%> + <c£’ + ) (7] + f 2 <T )] <fe* 

' 1 >A C = “X 1 [( <l'e + ) ( ■ f[ + <t>r 2 )5‘! 


j(t) 

J 33 


* =1 h_ 


~ h (a™C™ + a£>C<? + ) (t; + 7^) 


5 3 


<7z. 


(B15d) 


(B15e) 


<1>A C = “Z I [(?«* +^"J( r i +^2^ 


(*) 

33 


* =1 Vi 


- * K’d? + + <c£> ) ( 7 ; + u) 


5 3 


dz. 


(B15f) 


mM l = “X J [(?7* + ^7v)( r i + M ) S \ 

k=1 K-y 


(*) 

33 


)(T x +T^ 

■)f£ 

<r 3 V 

/ V 1 

v 

3 J. 


tfe. 


(B15g) 


Appendix C 

The matrices b and P in Eq. (33) can be expressed as 


where 


V) a T = 

U C7 


) r (0) r (0) r (nr (i) r 

a , a , a , a , a , 

_py’ zz ’ x_p ’ xx’ XV ’ 

R r (0) R r (0) R r (0) R r (i) R r (i) R r (n R r (0) r (0) t \ 

H vt 1 Mr 1 Mzz> M a ’ r xx’ Mvy’ Mzv Mvz ’ Mvz j 

(Cla) 

. f«)„r (t) n T (£) T (t) T (£) T 1 

\ U '^(0)> “ftfOM ^7(2)’ “17(3) > ^7(4)/ 

(Clb) 

E-, tl 

J3 

O 

E-n ti 

_b 

0 

0 
E^ tj 

a 

o' 

11 

_b 

0 

(Clc) 

_r«)f>r (OijJ 1 (P R r (<) R r V)pT \ 

1 M=,( 0)> P#7(l)> P#7(2)> p^7(3)> M=7(4)q 

(Cld) 

(0) R r _|(0) R r (oinr (0) R r 1 

Mzz — \ Mzz(O) ’ Mrz(l)’ Mzz(2) j 

(Cle) 

h dT 

O 

O 

II 

^d J 

0 

(Clf) 


and 
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>N L o o 

0 0 


0 ""Q' 

) 0 0 

(') N r — / (, *lV r (') N r (/) N r ( 0 N r \ 

^f!7(0)» ^^7(1)’ ^#7(2)’ ^7(3)’ ^#7(4)/ 

( 0 ) N r _f( 0 ) N r ( 0 ) N r ( 0 ) N r 1 

^zz - ) i, z;(0)’ zz(l)’ i, zz(2)j 


(,) M r = l (,, M r ‘'’M 7, (,, M r (,) M r (,, M r 

iV1 f7 l 0)’ iV1 ^(2)’ 1V1 f7( 3)’ 1V1 f7(4) 

<°>M r - / <0) M r (0) M r (0 ’M r \ 

lrl zz [ 1T1 -ZZ(0)’ iT1 -zz(l)> iT1 -zz(2)j 


(°)n r -!< 0 >n T (0) n T ( °>n T (0, n T 

Vf, \ V^( 0 )> V^(l)> V|^(2)> Vf,( 3 ) 


The matrix L' ( 1 in Eq. (37), which consists of the interpolation functions, is defined as 


~N,(s) N 2 (s) 0 0 


N,(s) N 2 (s) 0 


L' (i, = 


H t (s) H 2 (s) 


HJs) H Zs (s) 


0 0 
0 0 


HM) H , (.s) 




^,(■0 H 2 (s) 


3M 2 (s) 3 M 3 (s) 


6M } (s) 6M,(s) 


(continued on next page) 
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0 0 

0 0 

0 0 

0 0 

0 0 

N,(s) N 2 (s) 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 


0 

0 

H 3 (s) 

0 

0 

0 

H Xs (s) 

0 

0 

0 

0 

3 M 3 (s) 

4 

0 


0 0 0 

0 0 0 

H 4 (s) 0 0 

0 0 0 

0 0 0 

0 0 0 

H 4s (s) 0 0 

0 N t (s) N 2 (s ) 

0 0 

0 0 

0 0 

0 0 

0 0 


0 

0 

0 

3 M 3 (s) 

4 

0 


0 0 

0 0 

0 0 

H 3 (S) H 4 (s) 

0 0 

0 0 

0 0 

0 0 

# 3 ,» H 4 Js) 
0 0 

0 0 

0 0 

0 0 


0 0 

0 0 

0 0 

0 0 

H 3 (S) H 4 (s) 

0 0 

0 0 

0 0 

0 0 

0 0 

H 3 ,» H 4s { S ) 
3 M } (s) 3 M } (s) 

5 5 

0 0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

,, , , M 3 (s) 
0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

,,, , M 3 (s) 
0 


0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

N,(S) N 2 (s) 


0 

0 

0 

0 

0 

0 

0 

0 

0 

N,(s) 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

N 2 (s) 

0 

0 

0 


(C3) 


The matrix Z u) in Eq. (39), which is the transformation matrix between the nodal degrees of freedom defined in the 
element (x, y ) and line coordinate («, s) systems for the k th edge, is defined as 


[Zi,Z 2 , 0 ] 

(C4a) 

[0,z„z 2 ] 

(C4b) 

[z 2 ,o,z,] 

(C4c) 


where Z, and Z 2 are defined below in terms of c t and s h representing c k = cos((/9 ; ) and s k ~ sintr/?, ) , respectively 
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The transformation between the displacement vectors, u' w and u (t) , is established as u' (t ’ ) =Y (ir) u (t> , where the 
matrix Y U ) is defined as 

" c k s k 0 0 0 0 0 0 0 0 0 0 0 " 

-s k c k 0 0 0 0 0 0 0 0 0 0 0 

0010000000000 
0001000000000 
0000100000000 
0 0 0 0 0 c k s k 0 0 0 0 0 0 

Y (k) = 0 0 0 0 0 ~s k c k 0 0 0 0 0 0 (C5) 

0 00000 0 c, i t 0 0 0 0 

0 0 0 0 0 0 0 -s k c k 0 0 0 0 

0 00000 0 0 0 c k s k 0 0 

0 0 0 0 0 0 0 0 0 -s k c k 0 0 

0 0 0 0 0 0 0 0 0 0 0 c k s k 

0 0 0 0 0 0 0 0 0 0 0 -s k c k 
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The matrix B u in Eq. (47) is defined as 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

h 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

h 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

5/4 
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0 
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0 

0 

5/4 
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0 

0 

0 

0 

5/4 

0 

0 

0 

1 

5/4 

0 


(C6) 


The matrix B s in Eq. (48), which relates the boundary stress vector T' h l:> to the resultant stress vector defined in the 
element, is defined as 


B, = 
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0 
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0 


0 


0 

0 

0 

0 

0 

0 

n 

0 

", 
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0 

0 

0 

0 
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0 

0 
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", 
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0 

0 


0 


0 

", 

0 

0 

« 

0 

0 

0 

o -/?j 

' 3/ 
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v dx) 

-/? 

f d ) 

n y yr 
\ ' dy) 

-h\ 

f 3 d ) 

n — + /7 V — 

v ' dx dy j 

0 

0 

0 

0 

0 

0 

0 

0 

", 


0 


"v 
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0 

0 

0 

0 

0 

0 

0 

0 
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0 

0 

0 

0 

0 

0 

0 
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0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 


0 

0 

0 

0 

0 




0 

0 



0 

0 

o" 








0 

0 



0 
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r 3 / 
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d) 

5 ( 

d 

. 3 ^ 







"ir 

0 
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0 

0 

0 

-« v 


dx , 


dy 


dx 


-n 


«v 

0 


0 

0 

0 

0 

0 

( d_ 
' y dy 

0 

n„ 


+ ", 

dx dy 

0 

0 

0 

0 

0 


n — + n — | 0 0 

3 dx x dy 


n. 


(C7) 
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The components of the resultant stress vector, S h , given in Eq. (29) satisfy the equilibrium equations provided 
that the following constraint conditions are imposed in order to eliminate the redundancy of the assumed stress field: 


Re 


Re 


Y J z k -'(k i0) a xxm +ik (0) a xym ) 
.* = 1 


= 0 


= 0 


Re 


( 4 


X z*- 2 (k(k - 1) (1 7U, - *(* - !) °X(*) + 2ik ( k - D ) 


\k = 2 


( 3 

, N V 

Z-"( 

k m Pxzm +i k m P Mk )) J 


= 0 


Re 


; / 4 A . 

f 4 

A 

- U 


1 

-t 

S' 

V. 

U=1 

J 


i> a ( io vw 

\k=0 


= 0 


Re 


Z zA ( ( Vvz(t)) 

\k = 0 


= 0 


Re 


f 4 


X z*- 2 (*(* - 1) ^ - k(k - 1) + 2 ik(k - 1) <X <*) ) 

\k = 2 


1 

+ — 
h 


( 2 > 

k (0) 

z=(i) 

Vy=o y 


Z z a - 


= 0 


Re 


X^- 2 (^-i) (i x w -^-i) ( x (l , +2ik(k-ir^ m ) 

\k = 2 


1 


f 2 




iy (0, A.c*> 

V*=o y 


= 0 


(C8a) 

(C8b) 

(C8c) 


(C8d) 


(C8e) 


(C8f) 


(C8g) 


The relations used in the constmction of matrix c in Eq. (30) and also defined explicitly in terms of complex 
variables in Eq. (C8) are imposed as constraint equations that are suitable for computational implementation as 
follows: 


A:Re[ 

(0) a 1 

-A:Im[ 

- (0) a 

U xy(k) _ 

1 = » 

(k = 1,2, 3, 4) (C9a) 

A~Im^ 

m a 1 

^xi(k) J 

+ A~Re 

~ (0) a 

U xy(k) _ 

H 

(A: = 2,3,4) (C9b) 

-A:Im 

r (0) « 

L “xK*) 

] + Re[ 

~ m a 

U xy(k)_ 

H 

(k = 1,2, 3, 4) (C9c) 

A:Re[ 

’ (0, « 

U Mk)_ 

] + !m[ 

w a 1 

U xy(k) J 

= 0 

(A: = 2,3,4) (C9d) 
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(k = 2,3,4) (C9e) 


f (*(* - 1) Re[ (1, ^ w ] - k(k - 1) Ref ' ‘Xm ] - 2 k(k - 1) Im[<X ( J) 

-((*■ ~ 1) Re [ ( V ( *-.) ] ~ (k ~ 1) Im [ “Vv-a-o]) = 0 

f (*(* - 1) Im[ a, ^ (t) ]~k(k~ 1) !»[%,, ] + 2k(k - 1) Re [ ]) 

-((* - 1) H ' 1 V(*-d ] + (* - 1) Re [ (0 V^(*-i)]) = 0 

|(-A: Im [ (1 ^ { ^ ] + A: Re [ (1 >^ (t) ]) + (-A: Im [ ] + A: Re [ 

-Re[ , °V M ,- l) ] = 0 

^-f* Ref m Pyy( k ) ] + * Im[ ]) + (A: Re[ (0) ^ w ] + A: Im[ (0 Ak*) ]) 

- Im [ (0 V v . (t -_ 1) ] - 0 

|(^ Re[ ^ loaf ]) + (^ Re[ , ] - * I m [ (0> A : v W ]) 

-Re[ ( V, ( *-i)] = 0 

-H ro ^,)]=° 

\ (Re [ ( VJ) + h ( (*(* - 1) Re [ m a„ (t) ]-*(*- 1) Re [ m a )y , k) ] 
-2k(k - 1) Im [ (1) a^, (t) ] ) = 0 

^-(im [\ (H) ]) + *((*(* - 1) Im[ ( \ (t) ]- k(k - 1) Im [ \ (J) ] 
+2k(k - 1) Re [ (1 ) «*k*) ]) = 0 


(k = 3,4) (CA9f) 


(* = 1,2, 3, 4) (C9g) 


(* = 2,3,4) (C9h) 


(* = 1,2, 3, 4) (C9i) 


(* = 2,3,4) (C9j) 


(* = 2,3,4) (C9k) 


(* = 3,4) (C91) 


^r( Re [ (0) ^- 2) ]) + ((k(k - 1) Ref 1 'X*, ]-*(*- 1) Ref ( X t) ] 

-2*(*-l)Im[<X W ]) = 0 


(* = 2,3,4) (C9m) 


Y(H (0) ^ ( ,- 2 ) ]) + ((k(k - 1) Imf (1 X (t) ]-*(*- 1) Imf (1 X W ] 
+ 2*(*-l)Re[<% ( J) = 0 


(* = 3,4) (C9n) 


Appendix D 

The transverse shear stresses, (r[ k J and , in the * ,h layer are obtained by integrating the in-plane equilibrium 
equations while utilizing the constitutive relations given in Eq. (6) as 


<T™=- 


j( 

-A/2 


<T™+<T«\ 


)d£ 


(Dla) 
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°'-=- I «,+<’,)« 


yz J \ xy> x yy>y ) 

-h/2 


(Dlb) 


in which the derivatives of in-plane stresses are expressed as the dot products 


C7 a) = 

XX, X 

(c {k) 

r {k) 

^12 9 

R(k) 

7V 13 9 

c;r) 

■"y^xx,x 

-<XxAT x , 

£ >y,x 

~ a yy AT ,x> 

^ ‘ ZZ,X ’ 

Yxy,x 

~ a xy AT x ) 

(Die) 

a (k) =| 
yy,y ' 

(c {k) 

^12 9 

FKk) 
^22 ’ 

R (k) 

1V 23 9 

c«’)[ 


~ a xx AT ,y’ 

£ yy,y 

~ a yy AT ,y> 

® zz ,y 9 

Yxy.y 

~ a v AT .y) 

(Did) 

= 

xy,x 

(c (k) 

^16 9 

F*k) 

^26 9 

d(A) 

^63 9 

c“>) 

~( ^ xx,x 

-a xx AT x , 

£ yy,x 

~ a yy AT ,x’ 

^ ZZ,X 9 

Y xy,x 

Tx) 

(Die) 

a (k) =| 

xy,y \ 

\c {k) 
^16 ’ 

'‘"26 ’ 

iV 63 ’ 

cff) 

y^xx,y 

~ a xx AT ,y’ 

£ yy,y 

"■A' 1 ',- 

® zz,y 9 

Yxy.y 

~ a xy AT ,y) 

(Dlf) 


Note that the material properties C? and R] k) are functions of £/’ (thickness direction) only. Similarly, the 
transverse stress component, <j[ k J , is obtained by integrating the transverse equilibrium equation as 


<7 W =- f +a ik) 

^ I \ xz,x yz,y 

-h/2 


)d%= J J {<7^ +2 °v]xy+°w,yy)dVdZ 

-h/2 -h/2 


(D2a) 


in which the double derivatives of a/J , crj* , and cr/' are expressed as the dot products 


A k) = 

XX, XX 

(5! ”, 

r {k) 
^12 ’ 

R (k) 

JV 13 > 

CS’) 

(i'xx.xx 

~ a xx AT ,xx’ 

^ yy,xx 

~ a yy AT ,xx’ 

^ ZZ, XX 9 

Yxy,xx 

~ a xy AT ,xx) 

(D2b) 

A k) =| 

yy,yy ' 

(c<«. 

'“'22 9 

R(k) 

^ 23 9 

c !,*')[ 

If 

A 

- a xx AT ,yy> 

f 

yy,yy 

- a yy AT ,yy> 

® zz,yy 9 

Y xy,yy 

- a xy AT ,yy) 

(D2c) 

A k) = 

xy,xy 

(c™. 

^26 ’ 

iV 63 ’ 

c“>) 

“y^xx,^ 

~ a xx AT ,xy’ 

£ yy,xy 

~ a yy AT ,xy ’ 

^ZZ,Xf ’ 

y xy,xy 

“■■AT,.) 

(D2d) 


The integration of the in-plane and transverse equilibrium equations requires the continuous representation of the 
strain components and the transverse normal stress component so that their differentiation can be performed. The 
nodal degrees of freedom obtained from the solution of Eq. (61) are then utilized to approximate their continuous 
behavior in each element by employing linear interpolation functions, N m (x, y ) , as 


r„(.) „(«) ,,,(«) ...(e) ,.,(«) ...(e) ...(e) ...(e) ...(e) (e) (e) /3(e) a(e)l_Y V ) y 

l U > V ’ W > W 1 > W 2 ’ W ,x ’ W ,y ’ W l,x’ W l,)” W 2,x’ W 2,y’Vx ’ U y \ ~ N m A ■> > 

m = 1 

r,.(e) ,,(e) (e) (e) (e) (e) (e) (e) (e) (e) (e) /j(e) /)(e) 1 


(D3) 


in which the quantities with subscript (m) denote their value at node “m." 

However, for simply supported panels, their global behavior can be approximated in terms of Fourier series as 


« {x, y) = X X u ’ m ' cos ( a m x ) sin ( At) 

w=l n = l 

OO 00 

v(x,y) = Y J Y J V’ m sin (a m x) cos (ft„y) 

m- 1 n = 1 

oo oo 

[w(x, j), w, (x,y ) , w 2 (j£T,>-)] = ^ [if"'" , If/"" , W" m ] sin (a m x)sin (/?„ y) 

m = 1 /;=1 

oo oo 

[w, (x,j>),w 2>x (x, j)] = XX[^r>^r,^T]cos(a m x)sin(^„j;) 


(D4a) 

(D4b) 

(D4c) 

(D4d) 
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(D4e) 


r w 


■ y (x,y), w Uy (x,y), w 2 _ y (x, y)] = £ £ [w 7 , W£ , JT™ ] sin (a m x) cos (/?„y ) 


^ (x,y) = ZZ®'”" sin(« m x) c °s(/?„j;) (D4f) 

m - 1 «=1 

00 00 

0 y (x,y) = ZZ°r cos ( fl '» x ) sin (A)') ( D4 g) 

w=l «=1 

where or m and J3 n are defined as mx /(2a) and nx /(2b) , respectively. The unknown coefficients are calculated by 
employing Gauss integration techniques 


jj m " = 1 (x,. , y. ) cos ) sin [P n y t ) y/, 

Ob e =i i= i 

(D5a) 

1 N e G 

y„m _ — ^ Z A M v Ke) (x ( . , y i ) sin [a m x j ) cos (f) n y , ) y/, 
ab e= \ j=1 

(D5b) 

[w m ,wr,wr] = \ ii[w (e) {x,y)A e) (x,y)M e) { x ,y)] 

ob e=1 i= i 

sin («,„*,.) sin {P„y,)w, 

(D5c) 

[w”,w”,wz ] = (*> y)>< {*’ y)] 

ob e= i i= i 

cos(a m x l )sm(/3 n y l )i// l 

(D5d) 

[w™ ,w™ = ^rZZ[ w S e) { x >y)M% (x,y)M? y (x,y) 

ob e =i i=1 

si n ( «,„t) cos ( P„ y, ) ¥, 

1 (D5e) 

( x i > T, ) sin ( a m x i ) cos ( At, ) Vi 

ob e =j i= i 

(D5f) 

®r = “7 Z Z A(e) °T ( x > > T, ) cos ( a m x , ) sin ( P n y l ) yr i 
ob e=] i =1 

(D5g) 


where ;V is the number of elements, G is the number of Gauss points for each element, and t// is the weight of the 
numerical integration at the i th Gauss point. The area of e' h element, A (e] , is given as 


1 

1 

1 

.(e) 

'1 




>4 e) 



(D6) 


The Fourier coefficients are then used to obtain a continuous representation of the displacement components, using 
Eq. (D4), which can now be differentiated analytically. The second partial derivatives of the displacement 
components are used to calculate the first partial derivatives of the resultant strain and curvature components 
employing Eq. (5), which in turn are utilized to compute the transverse shear stress. Note that the computation of the 
transverse normal stress in Eq. (2a) requires the second derivative of the resultant strain and curvature components, 
therefore the third order derivatives of the displacement components in Eq. (D4) are evaluated. 

Through-the-thickness variations of the displacement components are computed from the strain displacement 
relation as 
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( a = x,y ) (Did) 


u 


(k) _ 
a 


S W 


-A/2 

z<*> 


1 W 


-u za )d^ + u a (x,y,-h/2) 

~ (*, JO - - w 2 a (I 2 - 1/5 ) y I + K a (x, j, -h/2) 


z (i > 

= f 4*^ + «- (*, T> - A/2 ) (D7b) 

-A/2 

where u a (x,y,-h/2) and u z (x,y,-h/2) are computed from the original displacement expression given in Eqs. (1) 
and (3). The strain components are obtained by the stress-strain relations. 
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